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We present a determinant quantum Monte Carlo study of the competition between instantaneous 
on-site Coulomb repulsion and retarded phonon-mediated attraction between electrons, as described 
by the two dimensional Hubbard-Holstein model. At half filling, we find a strong competition 
between antiferromagnetism (AFM) and charge density wave (CDW) order. We demonstrate that 
a simple picture of AFM-CDW competition that incorporates the phonon mediated attraction into 
an effective- U Hubbard model requires significant refinement. Specifically, retardation effects slow 
the onset of charge order, so that CDW order remains absent even when the effective U is negative. 
This delay opens a window where neither AFM nor CDW order is well established, and where there 
are signatures of a possible metallic phase. 
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The electron-phonon (el-ph) interaction is responsible 
for many phenomena in condensed matter physics, in- 
cluding charge density waves (CDWs) and conventional 
superconductivity. While the el-ph interaction is well un- 
derstood in metals, the role of phonons in strongly corre- 
lated systems is less clear, in part because the interplay of 
strong electron-electron (el-el) and el-ph interactions can 
lead to competing ordered phases. Despite its difficulty, 
this is an important problem to solve because multiple ex- 
perimental probes have detected signatures of significant 
lattice effects in strongly correlated materials. For ex- 
ample, in the cuprate high-temperature superconductors, 
angle-resolved photoemission spectroscopy (ARPES) has 
observed "kinks" in the band dispersion, which have been 
attributed to the el-ph interaction^ as well as small po- 
laron formation in undoped Ca2_ x Na ;E CuOCl2P^ Addi- 
tional evidence for a significant el-ph interaction include 
strong quasiparticle renormalizations detected by STM, 4 
and studies which have qualitatively reproduced optical 
conductivity peaks by including phonons!^ Besides the 
cuprates, other materials with both strong el-el and el-ph 
interactions include the manganites^ and fullerenesP 

On general grounds, two effects are expected when el- 
ph interactions are included in a system with strong el-el 
repulsion. The first is that the two interactions renormal- 
ize each other. The phonons mediate a retarded attrac- 
tive el-el interaction, thus reducing the effective Coloumb 
repulsion, while the el-el repulsion suppresses charge fluc- 
tuations, and hence the el-ph interaction, which couples 
to them. The second effect is a reduction in the quasipar- 
ticle weight due to additional scattering processes, which 
at large el-ph couplings can lead to a polaron crossover. 

A natural model for studying the interplay of the el- 



el and el-ph interactions is the Hubbard-Holstein (HH) 
model, which has been studied using various numerical 
approaches producing sometimes contradicting results. 
Within dynamical mean field theory (DMFT), the sup- 
pression of the el-ph interaction depends on the under- 
lying phase, and antiferromagnetic (AFM)-DMFT has 
found a moderate increase i n the critical el-ph coupling 
for small polaron formationPiSI i n contrast, diagram- 
matic quantum Monte Carlo work on the t-J-Holstein 
model found a reduction in the critical el-ph coupling 
needed for small polaron crossover.^ Dynamical clus- 
ter approximation (DCA) studies investigated the effect 
of phonons on the superconducting T c , and found that 
phonons suppress T c at small doping levelsp' however, 
including longer range hopping terms in the presence of 
phonons enhanced 

In addition to renormalization effects arising from the 
interplay of the el-el and el-ph interactions, competition 
between ordered phases can occur. On a two dimen- 
sional square lattice, at half filling the Hubbard and Hol- 
stein models have instabilities towards (7r/a, 7r/a) AFM 
and CDW orders, respectively; these phases compete in 
the HH model. Due to the many-body nature of the 
problem, compounded by the many degrees of freedom 
in the HH model, in general there is no exact solution. 
In one dimension, the HH phase diagram has been es- 
tablished via several numerical approaches, with an in- 
termediate metallic state between the AFM and CDW 
phasesP^ffil] The size of the metallic region grows with 
increasing phonon frequencyP^H^ A similar competition 
between AFM and CDW orders and phase diagram have 
been mapped out in infinite dimensions with DMFTP^U 
The AFM-CDW competition in two dimensions also has 
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been studied with perturbative^^H as we [\ as strong cou- 
pling 26 techniques. 

In this work we present a determinant quantum Monte 
Carlo (DQMC) study of the two dimensional single-band 
HH model at half filling. DQMC is a numerically ex- 
act method that treats el-el and el-ph interactions on an 
equal footing and non-perturbatively. A nonzero el-ph 
coupling introduces the fermion sign problem at half fill- 
ingEZl Nevertheless, simulations for all parameter ranges 
presented in this work can be done down to T = IT/40, 
where W is the non-interacting bandwidth. Significantly 
lower temperatures can be reached in some regimes. For 
details of the DQMC method, please refer to Refs. l2"5H5Ul 

The Hamiltonian for the single-band HH model is H = 
Hkin + Hi at + Hint, where 
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Here <...> denotes a sum over nearest neighbors, c i(7 
creates an electron with spin a at site i, hi„ — c ia c ia , 
t is the nearest neighbor hopping, is the phonon fre- 
quency, U is the el-el interaction strength, g is the el- 
ph interaction strength, and /i is the chemical potential, 
which is adjusted to maintain half filling. The dimen- 
sionless electron-phonon coupling constant is defined as 
A = g 2 /MQ 2 W. Throughout we take t = 1, M = 1, and 
a = 1 as our units of energy, mass, and length, respec- 
tively. 

We first study the spin and charge susceptibilities Xs 
and Xc which are given by 
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The susceptibilities at wavevector q = (tt, 7r) are shown 
in Fig. [T] for several values of U. With increasing A, \s 
decreases, signaling that the el-ph interaction reduces the 
strength of the effective el-el repulsion. This decrease in 
Xs occurs immediately with the inclusion of nonzero A 
for low to intermediate U, while for large U, suppres- 
sion of Xs does not occur until the el-ph coupling is fairly 
strong (A=0.5 for U=8t, and A > 1 for 17=10*). As Xs 
shrinks, Xc increases, indicating a clear competition be- 
tween the spin and charge orders. For all values of U 
considered here, Xc is negligible up to a Independent 
critical A, at which point it grows rapidly. However, for 
strong el-el interactions (U = 8t, 10*), Xc is still rela- 
tively small even at A = 1, due to the strong tendency 
toward AFM still present. Interestingly, rather than con- 
tinuously growing with A, the CDW susceptibility peaks 




FIG. 1. (a) Xs{^,^) an d (b) Xe{jf,Tr) for several U on an 
JV = 8 x 8 lattice. Inset of (b) shows Xs (dashed lines) and 
Xc (solid lines) at U — it for several lattice sizes. The error 
bars in the inset are suppressed for clarity. The remaining 
simulation parameters are: /3 = 4/t, At = 0.1/t, Q, = t. 
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FIG. 2. Structure factors (a) S 3 (7r,7r) and (b) S c (ir,n) for 
a U e s Hubbard model (black) and the U — St HH model 
with several phonon frequencies. The remaining simulation 
parameters are: N = 8 x 8, /3 = 4/t, At = 0.1/t. 



and then decreases, for U = 2t-6*. We attribute this be- 
havior to the finite CDW transition temperature in the 
HH model, which will be discussed in more detail below. 
The inset in Fig. [ljb) shows Xs and Xc for [7 = 4* for 
several lattice sizes, demonstrating that the lattice size 
has little effect on our conclusions. 

Since one of the effects of el-ph coupling is to reduce 
the effective strength of U, we investigate how well a 
U c e Hubbard model can describe the physics of the HH 
model. Integrating out the phonon field in the HH model 



3 



yields a dynamic el-el interaction: 



U eS (oj) = U 
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In the anti-adiabatic limit (SI— >oo) this interaction be- 
comes instantaneous, and reduces to the form U e g — 
U — XW. A frequency- independent U e g Hubbard model 
is often used to describe the HH model, even at finite 
SI. For example, DMFT studies have found that such 
an ap proac h captures the low-energy physics of the HH 
model.™ 

Fig. [2] compares the spin and charge structure factors 
Ss,c(q) =«\ c (q)dJ )C (q)> at q = (tt, tt) of a frequency- 
independent U c s Hubbard model and the U = 8t HH 
model at several phonon frequencies. Up to A « 0.25, 
S s (n,Tr) in the U e g and HH models agrees for all SI con- 
sidered. Beyond this point, S 8 (tt,tt) is suppressed more 
slowly in the HH model than in the U e g model, due to 
the retarded nature of the el-ph interaction captured in 
HH. As SI increases, the HH result comes closer to the 
C/ e ff result, until by SI — At, the two models agree within 
the error bars. We also considered other values of U (not 
shown), and found that for a given Si the difference be- 
tween the HH and U B g results grows as U increases. In 
contrast to S s (tt,tt), iS c (7T,7t) calculated in the U e g and 
HH models does not agree for any A. Rather, S c (ir, tt) im- 
mediately rises in the U e g model, while in the HH model 
it remains small until A ~ 0.75, and then rises sharply. 
As the phonon frequency increases, the HH and U e g re- 
sults get closer, although they are still inconsistent at 
SI = At. This result is generic; while the U e g Hubbard 
model has a CDW phase for any U e g that is negative, Xc 
remains suppressed well beyond the A at which U e g =0, 
as is clear in Fig.JTJb) where U e g —0 at A= 0.25, 0.5, and 
0.75 for U=2t, 4i, and 6i, respectively. 

An additional difference between the HH and U c g mod- 
els is the CDW transition temperature. In the HH model, 
while Tafm — in two dimensions due to the Mermin- 
Wagner theorem, Tcdw is finite because the order pa- 
rameter has two states. DQMC finite size scaling stud- 
ie d 31 ' 32 ^ the Holstein model found that t(3 C DW = 8-11 
for SI — t and A = 0.25. While we did not perform a scal- 
ing analysis for the HH model, we expect T C dw to be 
in the same temperature regime, because while the in- 
clusion of U in the HH model localizes carriers (which 
would lower Tcdw), it also pushes the CDW transition 
to a larger A (which would increase Tcdw)- In contrast, 
Tcdw = in the attractive-?/ Hubbard model. The 
sharply peaked nature of Xc in Fig. [IJb), differing from 
the slow evolution of x s , may be due to the proximity of 
the temperature t/3 = 4 to Tcdw- 

We now turn to the spectral properties of the HH 
model. To avoid analytic continuation, we focus on the 
spectral weight near the Fermi level, which is obtained 
from the imaginary time propagator via the relation^ 

(3C(k, r = /3/2) = | / du>A(k, u)g{ u , (3) (4) 
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FIG. 3. (a) Local phonon propagator /3D(r = 0, t = j3/2) 
for several values of U at /3 = 4/t. (b) Xc(^,n) for U = 2t 
and several values of /3. The remaining simulation parameters 
are: N = 8 x 8, Ar = 0.1/t, Q = t. 



where C and A are the propagator and spectral function, 
respectively, and g(oj,fi) = u>/ sinh(/3w/2) for bosons 
and = 1/ cosh(/?w/2) for fermions. At low tempera- 
ture, g(uj,f3) is sharply peaked about co = 0, so that 
A(k,uj = 0) dominates the integral. We consider the 
local propagator C(r = 0) = X^k^W' which is re- 
lated to the low energy projected density of states via 
N(0) PC(y = 0,t = /3/2)/tt. 



The phonon propagator, defined as Dij(r)= 
<T T Xi(T)Xj(0)>~<X> 2 , contains information on 
phonon softening at the CDW transition. In the Hol- 
stein model, the phonon spectral function is peaked at 
the bare phonon frequency ±51 in a system without el-ph 
coupling; el-ph interactions renormalize the phonon 
frequency and lead to spectral weight at other frequen- 
cies. In particular, the appearance of spectral weight at 
lu = indicates the development of a static CDW lattice 
distortion, which is revealed by j3D(r = 0, r = (3/2) 
(abbreviated as f3Dfj/ 2 ), as shown in Fig. [3| For low 
el-ph coupling, (3Dp/ 2 is negligible, since the system 
is far from the CDW state. It then increases at the 
same {/-dependent A at which Xc rapidly increases in 
Fig. [ljb). This phonon softening indicates that the 
CDW formation may have a Peierls-like origin, in which 
case the Fermi surface could be restored during the 
transition from an AFM to a CDW insulator. Fig. |3jb) 
shows Xc at U = 2t for several lower temperatures. 
With decreasing temperature, the rise in Xc sharpens 
dramatically and also shifts to lower A, appearing to 
asymptote towards a divergence in the susceptibility at 
low temperature around A = 0.3. We also note that the 
peak and subsequent decay in the CDW susceptibility 
discussed earlier appears robustly as a function of 
temperature. 

The electronic spectral weight near u) = also of- 
fers insight into the AFM-CDW transition. In this case, 
(3G(r = 0, t = f3/2) (abbreviated as (3Gp/ 2 ) distinguishes 
between insulating and metallic systems, being in the 
low temperature limit when a gap is present, and finite 
if a band disperses through the Fermi level. Fig. |4^a) 
shows f3Gp/ 2 for several values of U. For small A, f3Gp/ 2 
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FIG. 4. Local electron Green's function /3G(r = 0, r = /3/2) for (a) several values of U at /3 = 4/i, and (b) at U — 5t for 
several j3. (c) Average double occupancy <n^n^> at /3 = 4/t. The dashed line indicates where <n^n^>=0.25. (d) $ c — $ B at 
j3 = 4/i. The remaining simulation parameters are: iV = 8 x 8, At = 0.1/t, f2 = t. 



decreases with increasing U, indicating the opening of 
the Mott gap. As a function of increasing A, ftGp/2 falls 
for U=2t as the CDW gap develops. For U = At and 6t, 
initially grows as the el-ph interaction reduces the 
effective el-el repulsion and the Mott gap closes, and then 
decreases quickly as the CDW gap opens. For all these 
U, (3Gp/ 2 begins to fall at the same A at which Xc in- 
creases in Fig. [ljb) and the phonon softens in Fig. [3^a) . 
For U — 8t and lOt (3Gp/2 grows slowly with A as the 
Mott gap narrows. 

What can the peak in ftGp/2 at intermediate A in 
Fig. EJa) tell us about the AFM-CDW transition? In 
Fig. pffb), we plot (3Gp/ 2 at U = 5t for several temper- 
atures. As the temperature is lowered, /3G^/ 2 decreases 
at small and large A as the Mott and CDW gaps open, 
respectively. However, at intermediate el-ph coupling, 
PGp/2 actually grows with decreasing temperature, be- 
havior that could arise from an intervening metallic phase 
between the Mott and CDW insulating states. This in- 
crease in spectral weight at intermediate A was observed 
robustly for several lattice sizes and values of U. In addi- 
tion, the possible implication of a Fermi surface in the in- 
termediate state, from the phonon softening in Fig.ga), 
further supports the idea of an intermediate metallic 
phase. 

To further explore signatures of this possible metallic 
state, we plot the average double occupancy <n^n^> in 
Fig. gc). The double occupancy distinguishes between 



(%, 7t) AFM and CDW insulators, where it is and 0.5, 
respectively. In a metallic state (or an AFM-CDW coex- 
istence state), <n^n^>= 0.25. We find that <n^n^> in- 
creases smoothly with energy through 0.25, which is con- 
sistent with an intermediate metallic state, rather than 
a direct AFM-CDW transition at a critical A, where a 
sharp jump would be expected. While the transition from 
low to high double occupancy may sharpen as temper- 
ature is lowered, we note that in the range t/3 = 2 — 5, 
we found much less temperature dependence in <n^n^> 
than in other quantities considered in this paper. 

A finite temperature U — A phase diagram for f3 = 4/t, 
depicting the difference of the charge and spin order pa- 
rameters, $ c — $8, is shown in Fig. Kd). Here, the or- 
der parameters are defined as $ s = 2^j(^«t — f 1 ^) 2 /N 
and $ c = Y,icr(^i<r ~ l) 2 /N. Lines denoting U cS =0, 
<n^n^>= 0.25, and the peak in ftGp/ 2 are also included. 
The dominance of AFM and CDW orders at large U and 
A, respectively, is apparent. However, a sizable transition 
region, where $ c — $ s ~ is clearly visible. The lines U e s 
=0 and <n^n^>— 0.25 lie in the center of the transition 
region, while the peak in j3Gp/2 is toward the side dom- 
inated by spin order. The coincidence of multiple quan- 
tities consistent with a metallic state in the intermediate 
region of the phase diagram corroborate the case for the 
existence of such a phase. 

To summarize, in this work we demonstrated a strong 
competition between AFM and CDW phases in the two 
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dimensional single band HH model, and found evidence 
for a possible intermediate metallic regime existing be- 
tween the ordered phases. We investigated how well an 
effective-?/ Hubbard model can describe the physics of 
the HH model, and found that while in some regimes 
the two models give comparable results, in general the 
retarded nature of the el-ph interaction leads to signif- 
icant differences. The U-X phase diagram determined 
in our study is qualitatively similar to that found by 
low temperature numerical approaches, with the C/ e ff =0 
line dividing the regions of dominant spin and charge 
order parameters. We found evidence for an intermedi- 
ate metallic phase in two dimensions, similar to previous 
l-d results. E3HI2 The size of the intermediate metallic re- 
gion shrinks as the interaction strengths grow, which is 
consistent with Refs. H5HTT1 where the metallic phase is 



found to terminate at strong couplings. These findings 
contrast with the infinite dimensional DMFT results of 
Refs. and [231 where a direct order to order transition 
was found. Potential explorations for future work include 
studying the effect of phonon frequency on the interme- 
diate metallic state, and understanding more precisely 
where the metal-insulator transitions lie. 
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